function F = cal_firm_vars(param, const, eqm)

    % Useful function    
    firmvar = @(var) var(eqm.q_index)';
    
    
%% Quality, Productivity

    % Quality level
    firm.qlevel = firmvar(const.Qgrid);
    
    % Productivity
    firm.lnz = full(sum(const.lnzQ.*eqm.q_indicator, 2));

    
%% Export probabilities
    
    firm.prob_exp = full(sum(eqm.ExportProbQ.*eqm.q_indicator, 2));    
    firm.prob_nonexp = full(sum((1-eqm.ExportProbQ).*eqm.q_indicator, 2));
    
    
%% Export intensity    
    
    rv_rev = 1-eqm.D_F./eqm.D1;   
    firm.exp_intst_exp = firmvar(1-rv_rev);

        

%% Sales, Indegree, Outdegree

    % Log sales    
    firm.lnsales_exp = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi1));    
    firm.lnsales_nonexp = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi0)); 
        
%     % Number of suppliers
%     firm.indegree_exp = exp(firm.lnsales_exp).^(1/param.beta_m).*firmvar(const.C_m.*eqm.theta_m);
%     firm.indegree_nonexp = exp(firm.lnsales_nonexp).^(1/param.beta_m).*firmvar(const.C_m.*eqm.theta_m);
%        
%     % Number of customers
%     tempsum = sum(const.phi_v.*kron(ones(1,param.Q),eqm.theta_v'), 1);
%     firm.outdegree_exp = exp(firm.lnsales_exp).^(1/param.beta_v).*firmvar(eqm.rv_ads.*eqm.C_v1.*tempsum);
%     firm.outdegree_nonexp = exp(firm.lnsales_nonexp).^(1/param.beta_v).*firmvar(const.C_v0.*tempsum);
    
    
%% Wage, Employment

    % Average wage per worker
    firm.wage = firmvar(param.w.*const.wage_grid);    

    % Employment (number of workers)
    frac = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma;
    firm.employment_exp = frac.*exp(firm.lnsales_exp)./firm.wage;
    firm.employment_nonexp = frac.*exp(firm.lnsales_nonexp)./firm.wage;      
    
    
%% Average log wage of suppliers

%     % Unweighted average supplier wage
%     avg_unwgt_lnwageS_q = 1./eqm.V.*sum(const.phi_v.*kron(ones(param.Q,1), ...
%         (log(param.w)+const.lnwage_grid).*(const.C_v0.*eqm.Pi0.^(1/param.beta_v).*eqm.integ.E_zv0 ...
%         +eqm.rv_ads.*eqm.C_v1.*eqm.Pi1.^(1/param.beta_v).*eqm.integ.E_zv1)),2)';
%     firm.avg_unwgt_lnwageS = firmvar(avg_unwgt_lnwageS_q);
%     
%     % Expenditure weighted average supplier wage
%     avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c.^(param.sigma-1)./eqm.V.*sum(const.phi_y...
%         .*const.phi_v.*kron(ones(param.Q,1), (log(param.w)+const.lnwage_grid).*eqm.Psig),2)';    
%     firm.avg_wgt_lnwageS = firmvar(avg_wgt_lnwageS_q);
      
    
%% Fraction of higher quality match

%     % Fraction of higher (or same) quality supplier
%     flag_S_higher_q = zeros(param.Q, param.Q); %row is buyer, column is supplier
%     for i = 1:param.Q
%         flag_S_higher_q(i,:) = [zeros(1,i-1), ones(1,param.Q-i+1)];
%     end
%     frac_S_higher_q = sum(flag_S_higher_q.*const.phi_v.*repmat(eqm.Vbar,param.Q,1),2)'./eqm.V;
%     firm.frac_S_higher_q = firmvar(frac_S_higher_q);
%     
%     % Fraction of higher quality customer
%     flag_B_higher_q = zeros(param.Q, param.Q); %row is buyer, column is seller
%     for i = 1:param.Q
%         flag_B_higher_q(:,i) = [zeros(i-1,1); ones(param.Q-i+1,1)];
%     end    
%     tempmatrix = const.phi_v.*repmat(eqm.theta_v',1,param.Q);
%     frac_B_higher_q = sum(flag_B_higher_q.*tempmatrix,1)./sum(tempmatrix,1);        
%     firm.frac_B_higher_q = firmvar(frac_B_higher_q);
            

%% Return

    F = firm; 

    
end